################################
########TIES t-1 to t-10########
################################

#Figure A.4

#############	
#############	
### t – 1 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.02)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 1)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost+lag.lncapdist+lag.TIEScost:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc1.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 2 ###
#############
#############
############################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.035)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 2)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost2", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost2+lag.lncapdist+lag.TIEScost2:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc2.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost2", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 3 ###
#############
#############
############################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.03)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 3)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost3", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost3+lag.lncapdist+lag.TIEScost3:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc3.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost3", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 4 ###
#############
#############
############################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.02)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 4)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost4", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost4+lag.lncapdist+lag.TIEScost4:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc4.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost4", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 5 ###
#############
#############
############################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.02)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 5)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost5", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost5+lag.lncapdist+lag.TIEScost5:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc5.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost5", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 6 ###
#############
#############
############################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.02)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 6)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost6", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost6+lag.lncapdist+lag.TIEScost6:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc6.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost6", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 7 ###
#############
#############
############################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.02)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 7)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost7", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost7+lag.lncapdist+lag.TIEScost7:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc7.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost7", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 8 ###
#############
#############
############################
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.02)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 8)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost8", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost8+lag.lncapdist+lag.TIEScost8:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc8.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost8", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


#############
#############
### t – 9 ###
#############
#############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.05, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.02,.02)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 9)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost9", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost9+lag.lncapdist+lag.TIEScost9:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc9.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost9", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()


##############
##############
### t – 10 ###
##############
##############
felm_marginal_effects <- function(regression_model, data, treatment, moderator, treatment_translation, moderator_translation, dependent_variable_translation, alpha = 0.01, se = NULL) {
  library(ggplot2)
  library(ggthemes)
  library(gridExtra)
  
  ### defining function to get average marginal effects
  getmfx <- function(betas, data, treatment, moderator) {
    betas[treatment] + betas[paste0(treatment, ":", moderator)] * data[, moderator]
  }
  
  ### Defining function to analytically derive standard error for marginal effects
  getvarmfx <- function(my_vcov, data, treatment, moderator) {
    my_vcov[treatment, treatment] + data[, moderator]^2 * my_vcov[paste0(treatment, ":", moderator), paste0(treatment, ":", moderator)] + 2 * data[, moderator] * my_vcov[treatment, paste0(treatment, ":", moderator)]
  }
  
  ### constraining data to relevant variables
  data <- data[, c(treatment, moderator)]
  
  ### getting marginal effects
  data[, "marginal_effects"] <- getmfx(coef(regression_model), data, treatment, moderator)
  
  ### getting robust SEs
  if (is.null(se)) {
    data$se <- getvarmfx(regression_model$vcv, data, treatment, moderator)
  } else if (se == "clustered") {
    data$se <- getvarmfx(regression_model$clustervcv, data, treatment, moderator)
  } else if (se == "robust") {
    data$se <- getvarmfx(regression_model$robustvcv, data, treatment, moderator)
  }
  
  ### Getting CI bounds
  data[, "ci_lower"] <- data[, "marginal_effects"] - abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  data[, "ci_upper"] <- data[, "marginal_effects"] + abs(qt(alpha/2, regression_model$df, lower.tail = TRUE)) * sqrt(data$se)
  
  ### plotting marginal effects plot
  p_1 <- ggplot(data, aes_string(x = moderator)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), fill = "grey70", alpha = 0.4) +
    geom_line(aes(y = marginal_effects)) + ylim(-.01,.03)+
    theme_fivethirtyeight() +
    theme(plot.title = element_text(size = 11.5, hjust = 0.5), axis.title = element_text(size = 10)) +
    geom_rug() +
    xlab(moderator_translation) +
    ylab(paste("Marginal effect of",treatment_translation,"on",dependent_variable_translation)) +
    ggtitle("Average marginal effects (t - 10)")
  
  ### exporting plots as combined grob
  return(grid.arrange(p_1, ncol = 1))
}
#########################
###Subset NA Omit Data###
#########################
omit <- na.omit(NDs[,c("lag.lnpop_gpw_sum", "lnNTL","lag.lnbdist1",
                       "lag.lnexcluded",  "lag.lncapdist", "lag.lngas", "lag.lnoil",
                       "lag.lnNTL", "ccode", "year", "gid", "lag.p_polity2",
                       "lag.TIEScost10", "lag.lngcp")])
omit<-data.frame(omit)

#####################################
###Calculate for Sanctioned Groups###
#####################################
###Estimate model with NA omit data
felm4.1<- felm(lnNTL ~
                 lag.TIEScost10+lag.lncapdist+lag.TIEScost10:lag.lncapdist+lag.lngcp+
                 lag.lnNTL+lag.lnpop_gpw_sum+lag.lnexcluded+lag.lnbdist1+
                 lag.p_polity2+lag.lnoil+lag.lngas
               | 
                 factor(ccode)+factor(year)|0|gid, data=data.frame(omit))
summary(felm4.1)

###Create marginal effects plot
jpeg("TIESsanc10.jpeg")
felm_marginal_effects(regression_model = felm4.1, 
                      data = data.frame(omit), 
                      treatment = "lag.TIEScost10", moderator = "lag.lncapdist", 
                      se = "clustered",
                      treatment_translation = "Sanctions Cost", 
                      moderator_translation = "(ln)Capital Distance", 
                      dependent_variable_translation = "Subnational State Capacity")


dev.off()

